# Required packages ####
#install.packages("quantreg")
#install.packages("orthopolynom")
#install.packages("plotrix")
require(quantreg)
require(orthopolynom)

# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
CHwd<-file.path(directory, "Scripts/CH")
Datawd<-file.path(directory, "Processed")
Figwd<-file.path(directory, "Results/Figures")

# Functions #### 
#setwd(CHwd)
#source("CHfun.R")

polynomial.values <- function(polynomials, x, matrix=FALSE)
{
  ## Changed copy of polynomial.vales from orthopolynom in order
  ## to add matrix argument
  require(polynom)
  n <- length(polynomials)
  if(!matrix) {
    values <- vector(mode="list", length=n)
  } else {
    values <- matrix(ncol=n, nrow=length(x))
  }
  j <- 1
  while(j <= n) {
    if(!matrix) {
      values[[j]] <- predict(polynomials[[j]], x)
    } else {
      values[, j] <- predict(polynomials[[j]], x)
    }
    j <- j + 1
  }
  values
}

Legendre <- function(x, n, normalized=TRUE, intercept=FALSE)
{
  ## Create a design matrix for Legendre polynomials
  ## x - numeric
  ## n - see orthopolynom
  ## normalized - logical, see orthopolynom
  ## intercept - logical, add intercept
  tmp <- legendre.polynomials(n=n, normalized=normalized)
  if(!intercept) tmp <- tmp[2:length(tmp)]
  polynomial.values(polynomials=tmp, x=x, matrix=TRUE)
}

# Code begins ####
## Polynomial transformations from (-1,1) to (0,1) ####
n<-5
norma<-F
### Legendre ####
leg_fun<-polynomial.functions(legendre.polynomials(n=n, normalized=norma))
J <- 69
x <- -1+2*(1:J)/(J+1)

### Original base (-1,1)
{
  plot(x,leg_fun[[1]](x),type='l',ylim=c(-1,1),col='1')
  for (j in 2:(n+1)){
    c<-paste(j,sep="")
    lines(x,leg_fun[[j]](x),type='l',col=c)
  }
}

### Base on (0,1)
{
  p.list<-as.list(rep(0,length(leg_fun)))
  for (i in 1:length(leg_fun)){
    p.list[[i]]=as.function(alist(x = , leg_fun =  , leg_fun(2*x-1)))
  }
  
  taus <- 1:J/(J+1)
  
  plot(taus,p.list[[1]](taus,leg_fun[[1]]),type='l',col='1',ylim=c(-1,1),xlab="Quantile",ylab="")
  for (j in 2:(n+1)){
    c<-paste(j,sep="")
    lines(taus,p.list[[j]](taus,leg_fun[[j]]),type='l',col=c)
  }
  leg_name=''
  for (j in 1:(n+1)){
    pol<-paste("'p.list",j-1,"'",sep="")
    leg_name<-paste(leg_name,pol,sep=",")
  }
  leg_name<-substr(leg_name,2,nchar(leg_name))
  names <- eval(parse(text=paste('c(', leg_name, ')')))
  
  leg_col=''
  for (j in 1:n){
    color<-paste("'",j,"'",sep="")
    leg_col<-paste(leg_col,color,sep=",")
  }
  leg_col<-substr(leg_col,2,nchar(leg_col))
  colors <- eval(parse(text=paste('c(', leg_col, ')')))
  legend(0.6,-0.6, names, colors, ncol = 2, cex = 0.65)
}

### Matrix with polynomials
{
  X=''
  for (j in 1:(n+1)){
    v<-paste("p.list[[",j,"]](taus,leg_fun[[",j,"]])",sep="")
    X<-paste(X,v,sep=",")
  }
  X<-substr(X,2,nchar(X))
  X<-paste("cbind(",X,")")
  X<-eval(parse(text=X))
}

### Figure whit Legendre Base (0,1)
{
  line_type=''
  for (j in 1:(n+1)){
    line_type<-paste(line_type,j,sep=",")
  }
  line_type<-substr(line_type,2,nchar(line_type))
  line_type<-eval(parse(text=paste('c(', line_type, ')')))
  
  matplot(taus,X,type="l",ylim=c(-1,1),xlab="Quantile",ylab="")
  legend(0.6,-0.6, names,col=colors,lty=line_type, ncol = 2, cex = 0.65)
  
  setwd(Figwd)
  png("Legendre.png",width = 500, height = 500)
  matplot(taus,X,type="l",ylim=c(-1,1),xlab="Quantile",ylab="")
  legend(0.6,-0.6, names,col=colors,lty=line_type, ncol = 2, cex = 0.95)
  dev.off()
}

## Getting Quantile Regression Coefficients  ####
setwd(Datawd)
load("QR_mod.RData")

## Plot of approximations in Figure 1 ####
ii<-which(rownames(model0$coefficients)=="female")
setwd(Figwd)
for (n in c(2,6)) {
  png(paste('Approx_',n,'.png',sep = ""),width = 500, height = 500)
  par(mar=c(5.1, 4.1, 4.1, 2.1), mgp=c(3, 1, 0), las=0,oma=c(0.8,0.8,0.8,0.8))
  plot(smodel0, parm = ii, ols=F, xlab = "", col=c("darkgoldenrod1","gray88"),ylab = "", cex = 1, pch = 19,type="l",lwd = 2,main="")
  mtext("Quantile",side=1,line=1)
  options(scipen=10)
  yoff<-0.03
  xoff<-0.65
  es<-plotrix::emptyspace(taus,model0$coef[ii,])
  legend(es$x[1]-xoff,es$y[1]+yoff,
         legend=c(expression(hat(beta[j]) (tau)),"95% C.I.",expression(tilde(tilde(beta))[list(j,k)])), 
         bty = "n",
         col=c("darkgoldenrod1","gray88","black"), 
         lwd =c(2,8,2), 
         lty=c("solid","solid","dashed"), 
         x.intersp=c(1,1),
         ncol = 3, 
         cex = 1.1
  )
  
  Title=paste("Order K=",n)
  mtext(Title,side=3,line=-0.5,cex = 1.5)
  
  Mod<-lm(model0$coef[ii,] ~ Legendre(2*taus-1,n))
  y<-Mod$fit
  lines(taus,y,type="l",lty=2,col="black", lwd = 2)
  
  dev.off()
}